Applying to the LIGO data of gravitational wave, GW150914

last update: 9/17 (2018)


In [1]:
import DSP
using PyPlot

In [2]:
include("../juwvid.jl")


Out[2]:
juwvid

In [3]:
### Using LIGO data, provided by https://losc.ligo.org/events/GW150914/ 
data=readdlm("fig1-observed-H.txt"); # Hanford: Get from https://ligo.caltech.edu/WA
t=data[:,1]
y=data[:,2];
data=readdlm("fig1-observed-L.txt"); # Livingston: Get from https://ligo.caltech.edu/LA
t2=data[:,1]
y2=data[:,2];

In [4]:
fig=PyPlot.figure()
ax = fig[:add_subplot](1,2,1,aspect=0.05)
PyPlot.plot(t,y)
PyPlot.xlabel("time [s]")
PyPlot.xlim(0.25,0.46)
PyPlot.ylim(-1.5,1.5)
PyPlot.title("Hanford")
ax = fig[:add_subplot](1,2,2,aspect=0.05)
PyPlot.plot(t,y2)
PyPlot.xlim(0.25,0.46)
PyPlot.ylim(-1.5,1.5)
PyPlot.title("Livingston")
PyPlot.xlabel("time [s]")


Out[4]:
PyObject Text(0.5,24,'time [s]')

In [5]:
# compute frequency indices corresponding to 0-500 Hz 
dt=t[2]-t[1]
nufft=1024
juwutils.frequency_to_index([0.0,500.0], dt, length(y),nufft)


Out[5]:
2-element Array{Float64,1}:
   0.0  
 210.022

In [6]:
# define the frequency range
fin=collect(linspace(1.0,211.0,nufft));
# compute the analytic signals
z=DSP.Util.hilbert(y); 
z2=DSP.Util.hilbert(y2);

STFT


In [7]:
tfrst=stft.tfrstft(y);
tfrst2=stft.tfrstft(y2);


Use fft.
Use fft.

In [8]:
fig=PyPlot.figure()
ax = fig[:add_subplot](1,2,1)
a=juwplot.wtfrshow(abs.(tfrst).*abs.(tfrst),t[2]-t[1],t[1],t[end],NaN,NaN,11.0)
PyPlot.xlabel("time [s]")
PyPlot.ylabel("frequency [Hz]")
PyPlot.title("Hanford")
PyPlot.ylim(0,500)
ax = fig[:add_subplot](1,2,2)
a=juwplot.wtfrshow(abs.(tfrst2).*abs.(tfrst2),t[2]-t[1],t[1],t[end],NaN,NaN,11.0)
PyPlot.title("Livingston")
PyPlot.xlabel("time [s]")
PyPlot.ylim(0,500)


Out[8]:
(0, 500)

Wigner Ville Distribution


In [9]:
include("../cohenclass.jl")
tfrwv=cohenclass.tfrwv(z,NaN,NaN,fin,NaN,0,"nufft");
tfrwv2=cohenclass.tfrwv(z2,NaN,NaN,fin,NaN,0,"nufft");


Single Wigner Ville
WARNING: replacing module cohenclass
Use nufft.
Single Wigner Ville
Use nufft.

In [10]:
fig=PyPlot.figure()
ax = fig[:add_subplot](1,2,1)
a=juwplot.tfrshow(tfrwv,dt,t[1],t[end],fin[1],fin[end],0.7)
PyPlot.xlabel("time [s]")
PyPlot.ylabel("frequency [Hz]")
PyPlot.title("Hanford")
ax = fig[:add_subplot](1,2,2)
a=juwplot.tfrshow(tfrwv2,dt,t[1],t[end],fin[1],fin[end],0.7)
PyPlot.title("Livingston")
PyPlot.xlabel("time [s]")


Out[10]:
PyObject Text(0.5,24,'time [s]')

Pseudo Wigner Ville Distribution


In [11]:
tfrpwv=cohenclass.tfrpwv(z,NaN,NaN,fin,NaN,NaN,0,"nufft");
tfrpwv2=cohenclass.tfrpwv(z2,NaN,NaN,fin,NaN,NaN,0,"nufft");


Single pseudo Wigner Ville
Use nufft.
Single pseudo Wigner Ville
Use nufft.

In [12]:
fig=PyPlot.figure()
ax = fig[:add_subplot](1,2,1)
a=juwplot.tfrshow(tfrpwv,dt,t[1],t[end],fin[1],fin[end],0.7)
PyPlot.xlabel("time [s]")
PyPlot.ylabel("frequency [Hz]")
PyPlot.title("Hanford")
ax = fig[:add_subplot](1,2,2)
a=juwplot.tfrshow(tfrpwv2,dt,t[1],t[end],fin[1],fin[end],0.7)
PyPlot.title("Livingston")
PyPlot.xlabel("time [s]")


Out[12]:
PyObject Text(0.5,24,'time [s]')

L-wigner distribution


In [18]:
#alias free sm
afwv=smethod.tfrsm(y,NaN,4,NaN,2)
afwv2=smethod.tfrsm(y2,NaN,4,NaN,2)
#L wigner distribution (L=2)
tfrlw=lwigner.tfrlw2L(afwv,4);
tfrlw2=lwigner.tfrlw2L(afwv2,4);


Single S-method
Use fft.
Single S-method
Use fft.

In [19]:
fig=PyPlot.figure()
ax = fig[:add_subplot](1,2,1)
a=juwplot.wtfrshow(tfrlw,t[2]-t[1],t[1],t[end],NaN,NaN,11.0)
PyPlot.xlabel("time [s]")
PyPlot.ylabel("frequency [Hz]")
PyPlot.title("Hanford")
PyPlot.ylim(0,500)
ax = fig[:add_subplot](1,2,2)
a=juwplot.wtfrshow(tfrlw2,t[2]-t[1],t[1],t[end],NaN,NaN,11.0)
PyPlot.title("Livingston")
PyPlot.xlabel("time [s]")
PyPlot.ylim(0,500)


Out[19]:
(0, 500)

polynomial Wigner Ville distribution


In [20]:
tfrpo=polywv.tfrpowv(y,NaN,NaN,fin,8,8);
tfrpo2=polywv.tfrpowv(y2,NaN,NaN,fin,8,8);


Single S-method
Use nufft.
Single S-method
Use nufft.

In [21]:
fig=PyPlot.figure(figsize=(10,5))
ax = fig[:add_subplot](1,2,1)
a=juwplot.wtfrshow(log10.(complex(tfrpo)),dt,t[1],t[end],fin[1],fin[end],0.7,"CMRmap",4,10)
PyPlot.colorbar(a,shrink=0.3)
PyPlot.xlabel("time [s]")
PyPlot.ylabel("frequency [Hz]")
PyPlot.title("Hanford")
ax = fig[:add_subplot](1,2,2)
a=juwplot.wtfrshow(log10.(complex(tfrpo2)),dt,t[1],t[end],fin[1],fin[end],0.7,"CMRmap",4,10)
PyPlot.colorbar(a,shrink=0.3)
PyPlot.title("Livingston")
PyPlot.xlabel("time [s]")


Out[21]:
PyObject Text(0.5,24,'time [s]')

In [22]:
fig=PyPlot.figure()
ax = fig[:add_subplot](1,2,1)
a=juwplot.wtfrshow(tfrpo./maximum(tfrpo,1),dt,t[1],t[end],fin[1],fin[end],0.7)
PyPlot.xlabel("time [s]")
PyPlot.ylabel("frequency [Hz]")
PyPlot.title("Hanford")
ax = fig[:add_subplot](1,2,2)
a=juwplot.wtfrshow(tfrpo2./maximum(tfrpo2,1),dt,t[1],t[end],fin[1],fin[end],0.7)
PyPlot.title("Livingston")
PyPlot.xlabel("time [s]")


Out[22]:
PyObject Text(0.5,24,'time [s]')

In [ ]: